/*
 *	Equation solver
 *
 *	Copyright (C) 2011 XiaoJSoft Studio. All Rights Reserved.
 *	Copyright (C) Ji WenCong <whs_jwc@163.com>
 *
 *	This program is free software: you can redistribute it and/or modify
 *	it under the terms of the GNU General Public License as published by
 *	the Free Software Foundation, either version 3 of the License, or
 *	(at your option) any later version.
 *
 *	This program is distributed in the hope that it will be useful,
 *	but WITHOUT ANY WARRANTY; without even the implied warranty of
 *	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *	GNU General Public License for more details.
 *
 *	You should have received a copy of the GNU General Public License
 *	along with this program.  If not, see <http://www.gnu.org/licenses/>.
 *
 */

#include <stdio.h>
#include <memory.h>
#include <malloc.h>
#include "../include/equation.h"
#include "../include/mmdiv.h"

/*
 *	fact* solveEquations(fact *equs, fact *consts, int maxX, int maxY);
 *
 *	Returns values of a group of equations.
 *
 *	*equs	->	the equations,
 *	*consts	->	the constants,
 *	maxX		->	the number of unknowns,
 *	maxY		->	the number of equations.
 */
fact* solveEquations(fact *equs, fact *consts, int maxX, int maxY)
{
	fact cstCounter, equFirst;
	int bExch;
	fact *callback;
	fact *comEqu, *comCst;
	fact frDivSet;
	fact frExch;
	int dx, dy;
	fact *ret;

	/*  Exchange equations  */
	bExch = 0;
	for (dy = 1; dy <= maxY; dy++) {
		if ((equs + getPositionXY(1, dy, maxX) - 1)->numerator != 0) {
			if (dy == 1) {
				goto Exch;
			}
			for (dx = 1; dx <= maxX; dx++) {
				frExch = *(equs + getPositionXY(dx, dy, maxX) - 1);
				*(equs + getPositionXY(dx, dy, maxX) - 1) = *(equs + dx - 1);
				*(equs + dx - 1) = frExch;
			}
			frExch = *consts;
			*consts = *(consts + dy - 1);
			*(consts + dy - 1) = frExch;
Exch:
			bExch = 1;
			break;
		}
	}

	/*  Nothing to exchange, exit...  */
	if (bExch == 0) {
		return(NULL);
	}
	if (maxX == 1) {
		/*  Only one unknown, solving  */

		/*  Allocating enough memory areas to contain results  */
		ret = (fact*) malloc(1 * sizeof(fact));
		if (ret == NULL) {
			return(NULL);
		}
		*ret = fractionDivision(*consts, *equs);
		for (dy = 2; dy <= maxY; dy++) {
			if ((equs + dy - 1)->numerator == 0) {
				continue;
			}
			if (fractionCompare(fractionDivision(*(consts + dy - 1), *(equs + dy - 1)), *ret) != 0) {
				free(ret);
				return(NULL);
			}
		}
	} else {
		/*  One more unknown, eliminating  */

		/*  Allocating enough memory areas to contain processing equations  */
		comEqu = (fact*) malloc((maxX - 1) * (maxY - 1) * sizeof(fact));
		comCst = (fact*) malloc((maxY - 1) * sizeof(fact));
		if (comEqu == NULL || comCst == NULL) {
			return(NULL);
		}

		/*  Filling memory areas with zero  */
		memset(comEqu, 0, (maxX - 1) * (maxY - 1) * sizeof(fact));
		memset(comCst, 0, (maxY - 1) * sizeof(fact));

		/*  Building simplified equation  */
		for (dy = 1; dy <= maxY; dy++) {
			if ((equs + getPositionXY(1, dy, maxX) - 1)->numerator == 0) {
				/*  Need not to eliminate, exit...  */
				continue;
			}

			/*  Build the fraction which we need to simplify equations  */
			frDivSet = fractionDivision(fractionCreate(1, 1), *(equs + getPositionXY(1, dy, maxX) - 1));

			/*  Simplifying equations  */
			for (dx = 1; dx <= maxX; dx++) {
				*(equs + getPositionXY(dx, dy, maxX) - 1) = fractionMultiplination(*(equs + getPositionXY(dx, dy, maxX) - 1), frDivSet);
			}
			*(consts + dy - 1) = fractionMultiplination(*(consts + dy - 1), frDivSet);
		}

		/*  Eliminating  */
		for (dy = 2; dy <= maxY; dy++) {
			if ((equs + getPositionXY(1, dy, maxX) - 1)->numerator == 0) {
				for (dx = 2; dx <= maxX; dx++) {
					*(comEqu + getPositionXY(dx - 1, dy - 1, maxX - 1) - 1) = *(equs + getPositionXY(dx, dy, maxX) - 1);
				}
				*(comCst + dy - 2) = *(consts + dy - 1);
			} else {
				for (dx = 2; dx <= maxX; dx++) {
					*(comEqu + getPositionXY(dx - 1, dy - 1, maxX - 1) - 1) = fractionMinus(*(equs + dx - 1),  *(equs + getPositionXY(dx, dy, maxX) - 1));
				}
				*(comCst + dy - 2) = fractionMinus(*consts, *(consts + dy - 1));
			}
		}

		/*  Get the results of new equations  */
		callback = solveEquations(comEqu, comCst, maxX - 1, maxY - 1);
		if (callback == NULL) {
			free(comEqu);
			free(comCst);
			return(NULL);
		}

		/*  Get the result of the first unknown in equations  */
		cstCounter = *consts;
		for (dx = 2; dx <= maxX; dx++) {
			cstCounter = fractionMinus(cstCounter, fractionMultiplination(*(equs + dx - 1), *(callback + dx - 2)));
		}
		equFirst = fractionDivision(cstCounter, *equs);

		/*  Allocating enough memory areas to contain results  */
		ret = (fact*) malloc(maxX * sizeof(fact));
		if (ret == NULL) {
			free(callback);
			free(comEqu);
			free(comCst);
			return(NULL);
		}

		/*  Setting results  */
		*ret = equFirst;
		for (dx = 2; dx <= maxX; dx++) {
			*(ret + dx - 1) = *(callback + dx - 2);
		}
		free(comEqu);
		free(comCst);
		free(callback);
	}

	return(ret);
}

/*
 *	int getPositionXY(int dx, int dy, int maxX);
 *
 *	Returns the specific location of a member of a two-dimensional
 *   array in a heap.
 *
 *	dx		->	X-axis coordinate of the member,
 *	dy		->	Y-axis coordinate of the member,
 *	maxX		->	Maximum length of x-axis.
 */
int getPositionXY(int dx, int dy, int maxX)
{
	return((dy - 1) * maxX + dx);
}

